% Matlab Code for replication of "Optimal Debt Maturity and Firm Investment"
% by J. Jungherr and I. Schott, September 2020
function [funval, output] = firm_optim(val,firmstate,PARMS,VAL,doall)
% Inputs:  Guess (K, Leverage)
% Take as given: the firm's state today

debt        = firmstate.DEBT;  
zz          = firmstate.zz;  

% Get parameters from PARMS
W           = PARMS.W;  
SDF         = PARMS.SDF; 
Z           = PARMS.Z;
lenz        = PARMS.lenz;
TAU         = PARMS.TAU;
DELTA       = PARMS.DELTA;
C           = PARMS.C;
PSI         = PARMS.PSI;
ZETA        = PARMS.ZETA;
GAMMA       = PARMS.GAMMA;
XI          = PARMS.XI; 
sigmaeps    = PARMS.sigmaeps;  
ETA         = PARMS.ETA;  
F           = PARMS.F;
rhoz        = PARMS.rhoz;  

if nargin <= 4
    doall           = 0;
end

% Assign inputs
kguess      = val(1); % guess for k
btildeguess = 0; % long-term debt = 0 in STD model!
bs          = val(1)*val(2); % gues for short-term debt
alpha       = PSI*ZETA;
beta        = (1-PSI)*ZETA;
Lstar       = (beta*Z(zz)*kguess^alpha./W).^(1/(1-beta));

futurez     = min(lenz,zz+1); % future Z is either lenz or zz + 1
cp          = 0; %
cv          = (1-PARMS.x)*(rhoz*VAL.vF(zz) + (1-rhoz)*VAL.vF(futurez));
Y           = Z(zz).*(kguess^PSI*Lstar.^(1-PSI)).^ZETA;

% Solve for EPSBAR
epsguess    = ((bs + GAMMA*btildeguess - kguess - cv)/(1-TAU) - (Y - F - W.*Lstar - DELTA*kguess - C*(bs + btildeguess)))/kguess;
% Define integral over epsilon
PHI_u       = 1/2 *(1+erf((epsguess-0)/(sigmaeps*sqrt(2)))); 
cond_mean_u = (0 - 1/(sqrt(2*pi))*exp(-1/2*(epsguess/sigmaeps).^2)) ./ (1/2 *(1+erf((epsguess/sigmaeps)/sqrt(2)))) * sigmaeps;
KTILDE      =  PHI_u*kguess + (1-TAU).*(Y - F - W.*Lstar - DELTA*kguess).*PHI_u + (1-TAU)*kguess*PHI_u.*cond_mean_u ;
scrapvalue  = (1-XI)*KTILDE;

% break-even condition for issuance for short-term debt
issuance_std = SDF.*((1-PHI_u)*bs*(1 + C) + bs*scrapvalue/(bs+btildeguess));

% break-even condition for issuance for long-term debt
issuance_ltd = SDF.*((1-PHI_u)*(btildeguess-debt).*(GAMMA + C + (1-GAMMA)*cp) + ...
    (btildeguess-debt)*scrapvalue/(bs+btildeguess));

% Costs of issuing debt
costs = ETA * (bs + max(0,btildeguess-debt))^2;

cond_mean   = (1/(sqrt(2*pi))*exp(-1/2*(epsguess/sigmaeps).^2))./(1 - (1/2 *(1+erf((epsguess/sigmaeps)/sqrt(2)))) )*sigmaeps;
QPRIME_val  = kguess*(1-TAU)*(1-PHI_u).*(cond_mean - epsguess);

% Compute long-term bond price 
P_ltd       = (SDF.*((1-PHI_u).*(GAMMA+C + (1-GAMMA)*cp) + scrapvalue/(bs+btildeguess)));
if btildeguess == 0
    P_ltd = PARMS.p_fl;
end

% Compute short-term bond price 
P_std       = (SDF .* ((1-PHI_u)*(1+C) + scrapvalue/(bs+btildeguess)));
if bs == 0
    P_std = 1;
end

% Compute etilde
eprime      = kguess - issuance_ltd - issuance_std + costs;

funval      = - (-eprime +                  (SDF.*QPRIME_val));

if eprime < 0
    funval =  99999;
end

if doall == 1
    % Define remaining outputs:
    output.P_ltd        = P_ltd;
    output.P_std        = P_std;
    output.epsbar       = epsguess;
    output.ktilde       = KTILDE;
    output.costs        = costs;  
    output.eprime       = eprime;
    output.Lstar        = Lstar;
    output.V            = -funval;   
    output.q            = QPRIME_val;
    output.Mstar        = 0;    
    output.Kstar        = kguess;
    output.Levstar      = val(2);
end
